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ABSTRACT 

An expansion of a density field or particle distribution in basis functions 
which solve the Poisson equation both provides an easily parallelized n-body 
force algorithm and simplifies perturbation theories. The expansion converges 
quickly and provides the highest computational advantage if the lowest-order 
potential-density pair in the basis looks like the unperturbed galaxy or stellar 
system. Unfortunately, there are only a handful of such basis in the literature 
which limits this advantage. This paper presents an algorithm for deriving 
these bases to match a wide variety of galaxy models. The method is based on 
efficient numerical solution of the Sturm-Liouville equation and can be used for 
any geometry with a separable Laplacian. 

Two cases are described in detail. First for the spherical case, the lowest 
order basis function pair may be chosen to be exactly that of the underlying 
model. The profile may be cuspy or have a core and truncated or of infinite 
extent. Secondly, the method yields a three-dimensional cylindrical basis 
appropriate for studying galaxian disks. In this case, the vertical and radial 
bases are coupled; the lowest order radial part of the basis function can be 
chosen to match the underlying profile only in the disk plane. Practically, this 
basis is still a very good match to the overall disk profile and converges in a 
small number of terms. 

Subject headings: methods: numerical — stellar dynamics — Galaxy: 
structure — galaxies: structure 



1. Introduction 

The basis function n-body force solver is optimal for studying the global 
response of galaxies to perturbations or stabihty (Earn & Sellwood 1995). This 
technique was developed for astrophysical problems by Glutton-Brock (1972, 
1973), Kalnajs (1976), Fridman & Polyachenko (1984) and more recently by 
Hernquist & Ostriker (1992) who dubbed it the self- consistent field (SGF) 
method. Orthogonal function expansions are attractive Poisson equation solvers 
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for two reasons: 1) the expansions can be chosen to filter the structure over an 
interesting range of scales and simultaneously suppress small-scale noise; and 
2) the algorithm is computationally efficient, scaling linearly with the number 
of particles. Mathematically, this entire class of algorithms relies on the general 
properties of the Sturm-Liouville equation (SLE) of which the Poisson equation 
is a particular case. This same approach is common in perturbation theories 
and so facilitates direct comparison between n-body simulation and linear 
perturbation theory. In addition, this approach is straightforward to parallelize 
(e.g. Hernquist, Sigurdsson & Bryan 1995); we find the algorithm scales linearly 
with the number of processors with low overhead. If the basis set resembles the 
equilibrium galaxy, most of the computational work is concentrated on resolving 
the perturbation rather than the equilibrium. 

This last point is also a disadvantage of this technique in applications 
to date. If the equilibrium does not look like the basis set, the technique 
becomes less efficient and noisy because the expansion series must be sufficiently 
long to represent the equilibrium even without the perturbation. This paper 
describes a general method based on a numerical construction of orthogonal 
bases which remedies this situation. Solutions to the fundamental equation, the 
Sturm-Liouville equation, are well-understood and well-behaved. A number of 
recently published algorithms take advantage of the special properties of this 
differential equation to yield high-accuracy solutions with low computational 
work. Harnessing these developments to our needs leads to an algorithm for 
computing orthogonal bases whose lowest-order function matches any given any 
regular equilibrium; spherical and three-dimensional cylindrical solutions are 
described in detail here. The basic algorithm will be described in §§. 

For the spherical case, the proposed algorithm is competitive in performance 
with evaluation by recursion relation used for the published bases cited above 
and has reproduced them with high accuracy as a check. The cylindrical basis is 
a bit more cumbersome: one may rely on the same numerical solution to tailor 
the basis in the radial or vertical direction but not both simultaneously. Here, 
I choose to derive the radial basis numerically. The lowest-order radial basis 
functions then take the form /(r) exp(±zfcz).Q These may then be adapted to 
the background by principal component analysis. Although more cumbersome 
to implement and more time consuming to execute than the spherical case, it is 
still fast relative to non-expansion-based solvers. The details of the cylindrical 



Bases resulting the other choice has been explored by Earn (1996) using a different approach. 
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basis are given in §|3.2.1. 



2. The algorithm 
2.1. Motivation 

Here, I will explicitly describe the spherical and three-dimensional disk cases 
but all others are analogously derived with little change. 

The Poisson equation separates in any conic coordinate system. Choice 
of separation constants gives a differential equation in the SLE form for each 
dimension. The simplest solution employs the eigenfunctions of the Laplacian 
directly. For example, consider an expansion in spherical polar coordinates. 
Assuming that the density is proportional to the potential, the solution to 
Poisson's equation takes the form of an eigenfunction of the Laplacian: 

^ ^ 2<MM _ (([+11 

dr^ r dr 

The well-known full solution is the product of spherical harmonics in 9 and 
and Bessel functions in r. For a finite-radius mass distribution with an inner 
core, the inner boundary condition is the usual dR/dr\Q = and the multipole 
expansion provides the outgoing boundary condition: 



dR{r) 



dr 



rt 
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where is the outer edge of the profile. Using these boundary conditions and 
the orthogonality relation of the Bessel functions leads to the following potential 
and density pair: 
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where is the n*'^ zero of J/_i/2 and rt is the outer edge of the profile (Fridman 
& Polyachenko 1984). The functions and d^f^ have the following inner 
product: 

/•no 

rfrr2plr(r)<7*(r) = -5„„,. (4) 







Properties of solutions to the SLE ensure that this expansion set is complete. 
Therefore given a density distribution, the gravitational potential and force can 
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be found directly by expansion. The set d„) are often called biorthogonal. 
A similar expansion obtains for cylindrical polar coordinates. 

This straightforward approach has flaws. Bessel functions do not look 
like galaxian profiles and therefore accuracy demands high-order expansions. 
The required number of functions increases for extended profiles since Bessel 
functions are only orthogonal over a finite domain. To get around this, one may 
map the radial coordinate from the semi-infinite real axis to a finite segment. 
Appropriate choice of this transformation leads to new sets of biorthogonal 
functions in both the spherical (Clutton-Brock 1973, Hernquist & Ostriker 1992) 
and two-dimensional (Clutton-Brock 1972, Kalnajs 1976) and three-dimensional 
(Earn 1996) cylindrical cases. This small number of choices results in a 
mismatch between the lowest order basis functions and equilibrium profile. A 
poor fit between the basis and the underlying density profile is a source of noise 
in the force field which leads to relaxation (cf. Weinberg 1997). This is the 
general situation unless one's galaxy fortuitously coincides with particular sets 
of orthogonal polynomials or functions analytically derived from exact solutions 
of the Poisson equation. 

The solution proposed here is a numerical solution of the SLE using recently 
developed and published techniques (Marietta & Pryce 1991, Pruess & Fulton 
1993, see Pryce 1993 for a review). This allows adaptive construction of an 
expansion basis which matches the underlying density profile exactly and 
thereby removes one of the major limitations of this approach. The details are 
described in the next two sections. 

Alternative solutions to the mismatch problem have been described by Allen, 
Palmer & Papaloizou (1990) and Saha (1993). Both of these methods in their 
general form rely on the orthogonalization of a covering but non-orthogonal 
basis. There are two advantages to the approach developed here. First, the 
background profile is represented in one basis function with potentially rapid 
convergence in the perturbation. The basis evaluation is easily incorporated 
into existing SCF codes. Second, the same biorthogonal series may be used in 
linear perturbation analyses (e.g. Kalnajs 1976, Fridman & Polyachenko 1984, 
Weinberg 1990) and coefficients directly compared with n-body simulation. This 
development was motivated for precisely this reason and will underlie future 
inquiry. 
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2.2. Reduction of the Poisson equation to Sturm-Liouville form 

We present the cyhndrical polar case here to be explicit but again the 
others are analogous. The Laplace equation separates into the following three 
equations for a potential of the form \l/(r) = R{r)Z{z)Q{6): 



-—r—Rir) 
r dr dr 



m 



R{r) 
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Z(z) + PZ(z) 



Q{e) + m'Q{e) 






0. 



(5) 



Following the authors cited in |2.1| , we can look for a solution to the Poisson 
equation whose potential and density have the form 

^{r,z,e) = '^o{r)u{r)Z{z)e{e) 

p{r,z,d) = p,{r)u{r)Z{z)Q{d). (6) 
The Poisson equation then takes the form 
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together with second two of equation (^) above, where A is an unknown constant. 
The general form of the SLE is usually quoted as: 



d / du 
dx \ dx 



+ q{x)u = Xw{x)u 



where p{x),w{x) > over the domain of interest, [a,b]. The eigenf unctions are 
orthogonal (see Courant & Hilbert 1953 for extensive discussion) and may be 
normalized: jj^ dxw{x)u'^ = 1. Equation (|^) is easily rewritten in this form and 
one finds: 



d 
dr 



^2/ .du(r) 



m 



r"^ o{r)u{r) = 
A7rG\r'^o{r)po{r)u{r) (9) 

where Vr denotes the radial part of the Laplacian operator. The unknown 
constant A is the eigenvalue. Comparing to the standard SLE form, we have 
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These coefficient functions now provide tlie input to tlie standard packaged SLE 
solvers eitlier in tabular or subroutine form. The orthogonality condition for 
this case is 

-AttG drr-^o{r)po{r)u{ry = -AttG drr^p = l. (13) 

In other words, equations (H) are potential-density pairs. It is convenient to define 
p = AnGp so that the biorthogonality relation becomes J drr '$r{f)ps{i^) = — 5rs- 
Analogous expressions obtain for the spherical polar case. This development 
does not require that \l/o and po solve the Poisson equation but they must obey 
the appropriate boundary conditions at the center and at the edge (which may 
be r = oo). If we choose \E'o and po to be a solution of the Poisson equation then 
the lowest eigenvalue is unity and the eigenfunction u{r) is a constant function. 



2.3. Numerical solution 

For our problem, the SLE is well-conditioned and generally stable. Solutions 
may be straightforwardly obtained by shooting methods and standard ODE 
packages. Here, I used the Pruess method (Pruess 1973) as implemented by 
Pruess & Fulton (1993) with excellent success. Rather than find an approximate 
solution to the exact differential equation in the usual way, this approach 
approximates the differential equation by a piecewise continuous function — a 
discrete grid — and finds an exact solution to the approximate problem. The 
grid may be successively refined to ensure convergence to the desired tolerance. 
Additional numerical analysis provides the optimal choice of grid over the 
domain (which, again, may be infinite). This choice of a non-uniform grid is the 
numerical analog to transformation of the infinite interval to a discrete segment 
which plays a defining role in Glutton-Brock's approach. 

The resulting numerical eigenfunctions must be tabulated for future use. By 
contrast, the orthogonal polynomial schemes yield explicit recursion relations 
and this lack is the only practical disadvantage to this approach. On the other 
hand, the numerical SLE approach gives us the fiexibility to specify \E'o and po 
arbitrarily. For example, we may use the density profile from a previous n-body 
simulation. 



3. Examples and comparisons 
3.1. Spherical solutions for galaxian halos and spheroids 
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3.1.1. Method 

The boundary conditions must be appropriate for the problem at hand. In 
the case of spherical symmetry, there is a boundary at r = and r = rt. The 
inner boundary condition may be the traditional \E'' = or that for a scale-free 
cusp. The outer boundary condition follows from the multipole expansion: 

— = (14) 

dr r 

We may have — > oo in which case equation ([T^) applies in this limit. Once the 
functions are tabulated, the force algorithm proceeds as usual for an SCF code. 
Given <l>o and po, equations (p!OD-(p^) define the eigenvalue problem for the SLE. 
For example, the Pruess & Fulton code returns the eigenfunctions u{r) and the 
potential-density pairs follow from equations (^). The basis functions can be 
periodically recomputed to adaptively fit an evolving distribution; we have not 



implemented this for the spherical case here but see §3.2.1 and Weinberg (1996). 



3.1.2. Examples 

To test the spherical implementation, I assigned and po to the Hernquist 
model (Hernquist 1990) and compared the SLE solution with the analytic 
recursion relations (Hernquist & Ostriker 1992) for radial order n < IQ and 
m < 2. Performance of the spherical algorithm is well-documented so a 
comparison of potential pairs suffices. For m = 0, the numerically determined 
functions differed from the results of the recursion relation by one part in 10^ 
near the center and one part in 10^ elsewhere. This difference is due to the 
extrapolation of the cusp at r = 0. Here, the boundary condition for the cuspy 
profile fixes the asymptotic value of ratio ^'^/'^o as r — > 0. For m > the 
differences are obtained to the specified tolerance (one part in 10^ for these 
tests). To recover the Glutton-Brock (1973) set, one assigns and po according 
to the Plummer law; in this case, differences between the SLE solution and 
recursion relations are obtained for all m to the desired tolerance. In all cases, 
the orthogonality relation remains accurate and the potential density pair is an 
accurate solution of the Poisson equation. 

The background galaxian profile need not have finite mass and may be 
cuspy. For example, a basis set tailored to the singular isothermal sphere only 
requires one to specify appropriate boundary conditions. Boundary conditions 
corresponding to a disturbance not felt by in the singular core and at large radii 
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Fig. 1. — Potential-density pairs for Z = m = labeled by order, n = 1, . . .4 (upper and 
lower panels, resp.) whose lowest order member (n = 1) is the singular isothermal sphere. 
The density eigenf unctions are multiplied by r^. 




are: 

J d^{r)/dr = 
\^(r) = 

and 

jd'^{r)/dr = 1 = 

[{l + l)^{r)/r + d-^{r)/dr = 1^0 



(15) 



(16) 



where \l/(r) = \E'o(r)M(r). These same boundary conditions apply to the r^^'^ 
profile above. The / = boundary conditions ensure that the potential-density 
pairs are asymptotic to the spherical background at small and large radii. The 
/ 7^ boundary condition at small radius is the standard zero potential that 
ensures a single valued function. At large radius, we choose the condition 
obtained for an outer multipole. The four lowest-order / = pairs are shown 
in Figure |l|. The density functions are multiplied by oc l/po and, again, the 



-1 -0.5 0.5 1 

log ri/* 



Fig. 2. — Potential-density pairs for I = m = labeled by order, n = 1, . . .4 (upper and 
lower panels, resp.) whose lowest order member (n = 1) is the spherical deprojection of 
the r^/'' surface brightness law with -Re// = 1. To better represent the cuspy density profile 
graphically, the density eigenfunctions are shown relative to the deprojected r^/^ law. 

lowest order relative density function is constant as expected. 

In addition, the background galaxian profile need not have an analytic form. 
For example, the spherically symmetric profile that results in the empirical r^/^ 
surface density law may be numerically deprojected, tabulated and used as input 
to the SLE routines described above. A few of the lowest order potential-density 
pairs are shown in Figure ^. The density functions are shown relative to the 
background density. Notice that the lowest order relative density function is 
constant as expected. 



3.2. Three-dimensional cylindrical solutions for disks 
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3.2.1. Method 

For the cylindrical case, there are boundary conditions at r = 0, r = 
and z = ±zt. Here the situation is a bit trickier: the general solution requires 
matching outgoing boundary conditions in two dimensions. However as ^ oo, 
the multipole expansion implies that equation (|T^ applies to lowest order in 
1/r with / replaced by m. This technical simplification is strong motivation for 
adopting the radial domain r G [0, oo) as is done here. Implicit in equations 

and d^) is a separation constant chosen to give oscillatory functions Z{z) 
appropriate for a region of non-zero density. The functions match the outgoing 
Laplace solution at the outer boundary. By choosing the outer boundary of 
the 'pill box' sufficiently large (e.g. greater than ten scale heights), we obtain 
boundary conditions appropriate for the isolated disk. The vertical biorthogonal 
functions are then the sines and cosines of the discrete Fourier transform but 
over a vertical domain with twice the height of interest. This ensures that that 
force from density images do not affect potential (cf. Eastwood & Brownrigg 
1979). 

Experimentation suggests that 2® = 64 wave numbers are sufficient to 
adequately resolve the vertical structure. Separating real and imaginary parts 
(or equivalently, sine and cosine terms), this demands 128 coefficients per radial 
basis function! Although this trigonometric basis does not look the underlying 
basis, we can find an orthogonal transformation which rotates the basis into one 
which look like the desired equilibrium. We do this by an empirical orthogonal 
function analysis which is equivalent to principal component analysis (see 
Weinberg 1996 for details). In short, let the vector \E'j = {pij} be the potential 
basis functions evaluated at the position of the i*^ particle. The symmetric 
matrix S^^ = jj^ Y^iLi PitiPiv measures the weight of the particle distribution on 
the original basis. By diagonalizing this matrix, we determine an orthogonal 
transformation to a new basis. The lowest order basis function — the one 
with the largest eigenvalue — best represents the underlying point distribution 
followed in eigenvalue ranking by next best, etc. The first few functions usually 
represent most of the weight and this allows us to reduce the 128 coefficients to 
between two and six. 

Since the SLE solution is a good match to the radial profiles, we only need 
the empirical transformation in the ^;-direction. As an example of these new 
functions. Figure |^ shows the first three two-dimensional orthogonal functions 
for the two lowest radial orders based on a Monte Carlo realization of the 
exponential disk with unit scale length and scale height 1/10 using 10^ particles. 



Fig. 3. — Six orthogonal potential and density pairs (left and right panels, resp.) labeled 
by vertical index j and radial index n. Azimuthal order is m = 0. Five contour levels are 
linearly spaced from from zero to the largest absolute peak value. Positive (negative) levels 
are shown as solid (dotted) lines. 

Following the symmetry of the equilibrium model, the adaptive algorithm 
creates the lowest order modes with even symmetry about the disk midplane. 
However, the four or five lowest-order functions represent enough of the odd 
component to follow the evolution (cf. Fig. 

To summarize, the algorithm for the n-body force calculation for the 
three-dimensional cylindrical basis is then as follows: 

1. Compute S^i, from particle distribution using the basis derived from 
equation (P) with Z{z) chosen as discussed above. 

2. Compute transformation to new basis by solving for the eigenvectors. 

3. Retain eigenvectors corresponding to the M largest eigenvalues. The 
value of M may either be predetermined or computed adaptively from the 
cumulative distribution of eigenvalues (see Weinberg 1996 for details). 

4. Tabulate the new orthogonal set and use this to evaluate force for some 
time-interval on order of a dynamical time for the problem of interest. 

5. Goto 1. 

The computational bottleneck in this procedure is the construction of S^i, 
and computing Steps 1-3 can be a significant fraction of the integrated time 
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to advance the particles using the tabulated orthogonal functions for several 
hundred time steps (30% of the total for the case illustrated in Fig. ||). 
Nonetheless, the overall force evaluation is still very fast compared to other 
methods. 

Although the underlying trigonometric basis is bounded vertically from 
above and below, the boundary can be chosen large enough to permit arbitrarily 
large vertical distortions. Large vertical boundaries require more wavenumber to 
achieve a fixed resolution. In turn, more wavenumbers affect the computational 
overhead in computing the empirical basis but do not add to the CPU time 
required for the force evaluation itself. Therefore, large vertical boundaries 
remain practical as long as the transform to the empirical basis described in the 
algorithm above can be done infrequently. 



3.2.2. Examples 




Fig. 4. — First five density functions for m = (left) and m = 1 (right) with k = 27r/5 
ordered from bottom to top. The dotted curve on the lower-left-most panel shows the 
background exponential disk for comparison. 



Here we build a basis set that closely matches the typical exponential disk 
profile. As described in § p.2| , we adopt an axisymmetric separable density 
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coefficient amplitude 
0.5 




Fig. 5. — Expansion coefficient amplitudes for an exponential disk with sech^(z) vertical 
profile as a function of radial order and vertical wavenumber. 



profile, p{r,z) = pr{r)pz{z), chosen to match the background, Pr{r) = Poij). For 
this test, \E'o = — l/-y/l + {r/aY takes care of the boundary conditions. Recall 
that and po are not required to satisfy the Poisson equation; equation (|^ 
guarantees that the resulting basis functions will be orthogonal regardless. The 
results are shown in Figure ^ for the four lowest radial terms for m = and 1. 
The exponential scale length a = 1 and vertical boundary L = 10 is chosen to 
represent a disk with scale length to scale height ratio of 10:1. The wavenumbers 
are k = 27Tj/Lj,j = 0, 1, . . . ,jmax for piU box of half-height L. The density 
functions in the figure have = 27r/5 (j = 2). The lowest order m = case 
is compared with the exponential disk (dotted). For large k, the lowest order 
radial function falls off more rapidly than the exponential disk. However, the 
series converges quickly in radial order and wavenumber as demonstrated in 
Figure § which shows the coefficients for an expansion of Monte Carlo realization 
of an exponential disk with p^ oc sech.'^{z). Good agreement demonstrates that 
satisfactory results are obtained without exact Poisson solutions po and <l>o- The 
biorthogonality condition (eq. |T3p is good to one part in 10^. 
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The grid points for the Sturm-Liouville solution described in § p.3| are 
chosen by mapping the semi- infinite interval to the segment [—1,1] using 
X = (r — l)/(r + 1) and choosing points evenly spaced in x. The Pruess & 
Fulton algorithm can estimate the grid automatically to optimize accuracy but 
this mapping provided sufficiently high-accuracy and rapid execution. 

I checked accuracy and consistency of the final basis set by evaluating 
the gravitational force for a Monte Carlo distribution of 10^ bodies with the 
proposed method and with a direct summation. Contours of constant force are 
better than 1% except where the direct summation evaluation is badly affected 
by discreteness noise. 

4. Summary and conclusions 

This paper presents a numerical algorithm for constructing biorthogonal 
expansion bases for use in n-body force calculation and linear perturbation 
theory and explores its performance. The major results of this investigation are 
as follows: 

1. This algorithm removes one of the remaining limitations of the self- 
consistent ffeld (SCF) method by providing basis sets tailored to any 
background galaxian profile. 

2. The algorithm is applicable to any separable coordinate system. This 
paper details and benchmarks its implementation for spherical and 
three-dimensional cylindrical bases. 

3. Sturm-Liouville equation solvers are publically available (e.g. see Pruess & 
Fulton 1993 for Fortran code) and a desired basis is readily obtained using 
equations (|lO|)-(|ll). 

4. The main limitation of this method for n-body codes is the necessity to 
tabulate the basis functions rather than derive them from recursion relation 
on the fiy (as in Clutton-Brock 1973 and Hernquist & Ostriker 1992). On 
the other hand, this is largely a programming inconvenience; the algorithm 
is still easily parallelized and the table lookup is not a computational 
bottleneck. 

5. For spherical expansions, the algorithm is conceptually equivalent to and 
computationally competitive with the published SCF expansions. For 
three-dimensional cylindrical expansions, the coupling of the vertical 
and radial parts of the potential-density pairs requires an additional 
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orthogonalization step. This increases the computational overhead by up 
to 50% but does not effect scahng with particle number or parallelizability. 

6. The use of these basis sets is not limited to n-body simulation. They are 
easily used in semi- analytic linear perturbation calculations and, moreover, 
facilitate the comparison between the n-body and perturbation theory. 

I thank Lars Hcrnquiust and Neal Katz for discussion and suggestions. This 
work was supported in part by NSF grants AST-9529328 and the Alfred P. 
Sloan Foundation. 
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